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The surface scaling behavior of classical isotropic Heisenberg magnets is investigated by Monte - 
Carlo methods in d = 3 dimensions for various values of the surface - to - bulk coupling ratio Ji / J. 
For Ji/J < 1.0 critical behavior according to the ordinary surface universality class is found. New 
estimates for magnetic surface exponents are presented and compared to older estimates and their 
theoretical counterparts. For Ji/J > 2.0 scaling is still valid with effective exponents which depend 
on Ji/J. The surface structure factor S\(p,L) is investigated at bulk criticality as function of the 
momentum transfer p parallel to the surface and the system size L. For Ji/J < 1.0 and Ji/J > 2.0 
the full p dependence of Si (p, L) can be captured by generalized shape functions to a remarkable 
accuracy. Profiles of the magnetization and the energy density also confirm scaling, where for 
Ji/J < 1.0 the ordinary surface universality class is recovered and for Ji/J > 2.0 scaling with Ji/J 
dependent exponents is found. For Ji/J = 1.5 the system displays a striking crossover behavior 
from spurious long - range surface order to the ordinary surface universality class. For Ji/J > 2.0 
the effective scaling laws must be interpreted as nonasymptotic and the value Ji/J = 1.5 marks 
a crossover regime, in which the crossover from the nonasymptotic to the asymptotic (ordinary) 
surface scaling behavior can be resolved within numerically attainable system sizes. 
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I. INTRODUCTION 

Classical spin systems with O(N) symmetry display 
long - range magnetic order at low temperatures, if the 
spatial dimensionality is sufficiently high or the inter- 
action range is sufficiently large. Long - range order 
sets in spontaneously below a certain critical tempera- 
ture T — T c at zero magnetic field H — 0, where the 
field H has N components. In the vicinity of the critical 
point (T, H) = (T c , 0) these systems are characterized by 
an N component order parameter given by the magneti- 
zation. Both from the theoretical and the experimental 
point of view these systems provide important realiza- 
tions of the well known O(N) universality class for criti- 
cal behavior. In particular, the cases N = 1, N = 2, and 
N = 3 are most relevant for experiments, because they 
correspond to the Ising universality class (magnets with 
uniaxial anisotropies, simple and binary fluids), the XY 
universality class (magnets with planar anisotropies, 4 He 
near the normal - superfluid transition) , and the Heisen- 
berg universality class (isotropic magnets, e.g. Ni, EuO, 
EuS), respectively. Real samples of any material always 
have surfaces which display a distinct critical behavior 
in the vicinity of the critical point of the bulk material. 
In turn, the penetration depth of surface effects is set by 
the bulk correlation length £ = ^|r| _l/ which diverges 
at the critical point r = (T - T c ) /T c 0, where H = 
is assumed. The exponent v is the universal correlation 
length exponent and £± denotes the nonuniversal cor- 
relation length amplitudes above (+) and below (— ) T c 
(N = 1 only). It is now well established that surface 
critical behavior obeys the priciple of universality, i.e., 
critical exponents and scaling functions which character- 



ize surface critical behavior do not depend on microscopic 
details of the materialHtl General field - theoretic con- 
siderations based on the-i/) 4 Ginzburg - Landau model 
for semiinfinite systemsacl show that for O(N) symmet- 
ric systems there are three surface universality classes 
which are denoted as ordinary (O), surface - bulk or spe- 
cial (SB), and extraordinary (E). These three univer- 
sality classes can be briefly characterized by the absence 
or presence of surface order at the critical point of the 
bulk material. In particular, O stands for a disordered 
surface and E symbolizes an ordered surface, whereas 
SB indicates that surface and bulk are critical simultane- 
ously. The physical nature of the SB, surface universality 
class is that of a multicritical poinm. Surface order can 
arise spontaneously (E) or can be imposed from outside 
by a surface field Hi. The latter case is quite common 
in fluid systems with confining walls (Ising universality 
class, N = 1)q, where the bulk transition is then denoted 
as normal transition. It has been shown by rigorous argu- 
ments, that apart from corrections to scaling thejiormal 
and the extraordinary transition are equivalentti This 
is also quite important for the theory, because the nor- 
mal transition gives access to the extraordinary surface 
universality class even when the surface does not order 
spontaneously, i.e., the SB multicritical point does not 
exist. 

Within the framework of O(N) symmetric spin mod- 
els on lattices we will only consider the case of short - 
range interactions which is characterized by a nearest - 
neighbor coupling constant J in the bulk and by another 
nearest - neighbor coupling constant J\ in the surface. In 
this case the absence or presence of the SB transition cru- 
cially depends on the spatial dimension d and the value 
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of N. In d = 2 dimensions the SB transition does not 
exist, because the "surfaces" in this case are only one - 
dimensional. Note that for N > 2 a O(M) spin system 
does not show long - range order in d = 2u. In d = 3 the 
value of N becomes crucial. For N = 1 the SB transi- 
tion exists, because the surface, i.e., the Ising model in 
d = 2, exhibits long - range order below a finite critical 
temperature. In the phase diagram of a semi - infinite 
Ising ferromagnet the SB multicritical point is located at 
T = T C , H = Hi=0(7V=l here), and the critical value 
Jf/J of the surface - to - bulk coupling ratio Ji/J. Note 
that besides the lines of the O and E transitions a line 
of surface transitions terminates in the SB multicritical 
pointoH. For the Ising ferromagnet one finds Jf/J = 3/2 
on a simple cubic lattice to a remarkably high accuracy!! 
For N — 2 the surface of the spin system is a two - di- 
mensional XY model which does no longer show long - 
range orderQ. Instead, a line of Kosterlitz - Thouless sur- 
face transitions occurs in the (T, J\ / J) plane of the phase 
diagramH which terminates in a SB multicritical point at 
T = T c and J x j J = Jf/J ~ 1.5 for a XY ferromagnet 
on a simple cubic lattice. For N — 3 the SB transition 
does not exist, because the two - dimensional Heisenberg 
model has no phase transitionO. In the absence of symme- 
try breaking surface interactions the isotropic Heisenberg 
ferromagnet in d = 3 should therefore always display or- 
dinary surface critical behavior, regardless of the value 
of Ji/J. In d > 4 dimensions, i.e., within mean - field 
theory, the SB transition exists for any value of N. 

The dimensional crossover behavior of the SB multi- 
critical point outlined above greatly complicates the field 
- theoretic description of the SB transition for general 
N, apart from additional Goldstone problems below T c 
for N > lu which can be accountecUfpr in finite systems 
with periodic boundary conditional. Therefore, most 
of the field - theoretic work is devoted to the Ising uni- 
versality-class and to the ordinary transition for general 
N > lmLdO. In the following we will focus on the ordi- 
nary transition in isotropic Heisenberg magnets and we 
therefore only summarize some aspects of the ordinary 
surface universality class for later reference. For a recent 
general survey ftf surface critical behavior the reader is 
referred to Ref.cl. 

Apart from the standard bulk critical exponents addi- 
tional surface critical exponents must be introduced in 
order to describe critical behavior at the ordinary tran- 
sition. Let mi denote the modulus of the surface mag- 
netization per spin Mi /Area and let H and Hi denote 
the moduli of the bulk and surface field H and Hi, re- 
spectively. Then the critical behaviot, pif-Jpi and of the 
susceptibilities xi an d Xn is given bytrElO 
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for H = H± = and sufficiently small r. The sucepti- 
bilities \i an d Xn are known as the layer susceptibility 
and the surface or local susceptibility. In terms of the co- 
ordinates X|| parallel to the surface and z perpendicular 
to the surface ±he correlation decay exponents 771 1 and ?7_l 
are defined byaa 
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in the limit of large distances, where T = T c and H = 
Hi = is ass ume d. The surface exponents defined in 
Eqs.(l.l) and (F2)|-aije not independent. By virtue of 
the scaling relations 



/3i=i/(d-2 + T/||)/2 
71 = i/(2 - T] ± ) 

711 = 1^(1 - 7J||) 

v± = iv + v\\)/ 2 
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only one of the surface exponents defined above, say, r)\\ 
is independent for given v and 77, where 77 is the decay 
exponent of the bulk correlation function according to 

G (| x (l) - X ( 2 )|) ~ | x (l) _ x (2)|-(d-2+r ? ) at T= T ^ ^ 

remaining magnetic exponents 8\ and <$n defined by 
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are related to the bulk exponents v and r\ and the surface 
exponent t)\\ according to b~\ = v(d + 2 — rj)/(2(3\) and 
<5ii = v(d — 77||)/(2/3i), respectivelytlra. In analogy with 
Eq.(l.l) the layer specific heat C\ and the surface,or local 



specific heat C\\ display the critical singularitesa 

where e± is the surface energy density per spin and 
c = (Jf — Ji)/J is the surface enhancement. Within the 
framework of the field - theoretic renormalization group 
it has been shown by an explicit calculation of the energy 
density profildl3 and from the short - distance expansion 
of the energy density operator near the surfaceoEJ that 
Q'i = a — 1 and an = a — 2 — v at the ordinary transition. 
The critical singularities of all bulk and surface quanti- 
ties can therefore be expressed by the bulk exponents v 
and 77 and the surface exponent 7711, i.e., there is only 
a single independent surface exponent at the ordinary 
transitioned Within the field - theoretic renormalization 
group in d — A—e all surface exponents are known to two - 
loop order for general iVuO. The extrapolation of the e - 
expansion to d = 3 (e = 1) already yields reasonable esti- 
mates for the surface exponents, however, their numerical 
accuracy remains onknown due to the lack of higher or- 
der results which are needed to apply resummation tech- 
niques. Therefore, alternative analytic approaches have 
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been pursued in order to obtain improved estimates. In 
particular for thepHeisenberg model (3\ = 0.81 ± 0.04 
has been obtainedll3 which agrees with the extrapolated 
value for f3\U. The result of a perturbative calculation 
for the nonlinear sigma model in d = 2 + e has been 
combined with thc-results in d = 4 — e by means of a 
Pade approximanttll resulting in the numerical estimates 
/3i = 0.84 ± 0.01 and n\\ = 1.39 ± 0.02. A massive field- 
theoretic approach recently provided the new estimates 
Px = 0.880, 0.862, 0.889 from a [2/0], [0/2], and [1/1] 
Pade approximant in d = 3, respectivelyL3. From an 
earlier experiment on Ni(100) and Ni(110) surfaceslla the 
average estimates (3\ ~ 0.81 and fix ~ 0.79 have been 
obtained, respectively. 

For the. experimental verification of surface critical 
behaviora the strutuxe factor in s urfa ces and thin films is 
of particular interesttj. From Eq. ( L2 ) one concludes that 
for small parallel momentum transfer p the surface corre- 
lation function displays the algebraic p~ 1+ri w singularity, 
which gives access to the exponent rj\\ in a surface scat- 
tering experiment. Consequently, substantial theoretical 
effort has been spent on the theoretical understanding 
of correlation functions in surfaces and filmsE3 and the 
crossover behavior between themEd. A thorough survey 
of the properties of the static structure factor can be 
found in Ref.E3. Scattering experiments also give access 
to the order parameter profile which is governed by uni- 
versal shape functions. In d = 2 and at T = T c the shape 
functions of the order parameter and the energy densitiy 
profile for Ising and Potts models caa be obtained from 
conformal invariance considerations^. Generalizations 
to O(N) symmetric models in d = 2 are also possibles. 
In higher dimensions one has to resort to field - theo- 
retic methodsEj'Ej, where only the eseray density profile 
is nonzero at the ordinary transitiorli3€3. The identifica- 
tion of the surface universality class for a specific system 
is a rather delicate problem. For example, the presence 
of weak surface fields leads to a crossover from ordinary 
to extraordinary-behavior as the surface is approached 
from the interiorEH. For multi - component systems, e.g., 
binary alloys the surface universality classes ma;y-even be 
different for different orientations of the surfaceEa. 

Whereas most of the theoretical results quoted above 
can be applied to the case d = 3, N = 3 under consid- 
eration, numerical investigations have primarily focussed 
on the Ising universality class (N = 1) in d = 3. From a 
Monte - Carlo simulation of the isotropic Heisenberg fer- 
romagnet on a simple cubic lattice with an open surface 
on one side and a self - consistently determined surface 
field on the opposite sideEa an estimate fix = 0.75 ± 0.10 
has been found for the surface exponent fix (ordinary 
transition)Ej. The field dependenc e of the surface mag- 
netization nil at T = T c (see Eq.( |l.4| )) was determined 
lateitil and the estimate 8i — 2.3 ±0.1 was found. Trans- 
fer matrix Monte - Carlo calculations for the Heisen- 
berg ferromagnet,on small lattices provided the estimate 
fix = 0.80±0.03l3. Most of the numerical effort has been 



spent on the surface critical behavior of the Ising modelE 



and a few detailed studies also exist for the XY modelE-l 
The most accurate numerical estimates of Ising surface 
exponents and amplitude ratios at -the ordinary and the 
SB transition can be found in Ref.Q. We close this brief 
overview with the remark that the surface exponent fix 
for the ordinary transition is not affected by surface bond 
disorder or the presence of steps on the surfacea. A the- 
oretical explanation for this behavior can be found-from 
the construction of upper and lower bounds on 

Although a substantial wealth of information is already 
available for the surface critrical behavior of the isotropic 
ferromagnet with 0(N — 3) symmetry at the surface, the 
resulting estimates for the surface critical exponents are 
too disparate to provide a reliable basis for further inves- 
tigations, e.g, the surface contribution to the dynamic 
structure factor. Furthermore, a systematic scaling anal- 
ysis of static surface correlations is still missing, which 
provides essential information for the numerical analy- 
sis of experimental scattering data and Monte - Carlo 
data of the dynamic structure factor. Moreover, previ- 
ous investigations were focussed on the asymptotic scal- 
ing regime of the ordinary surface universality class. In 
real systems or computer simulations of isotropic Heisen- 
berg magnets, however, the surface properties at hand 
may be such that asymptotic scaling is only obtained 
after a crossover regime of a a certain width has been 
traversed. This also means that some sort of nonasymp- 
totic surface behavior must occur before the crossover 
regime is reached even if the bulk is already critical, i.e., 
displays asymptotic scaling behavior governed by bulk 
exponents. It is important for the analysis of numeri- 
cal data to localize these regimes and to describe their 
properties. 

It is the purpose of this investigation to fill at least 
some of the aforementioned gaps. In particular new in- 
dependent estimates for the surface critical exponents 
of the isotropic Heisenberg ferromaanet at the ordinary 
transition are provided. With RefsEa andc3 as guidelines 
for the universal form of the shape functions the surface 
structure factor and the order parameter and energy den- 
sity profiles are analyzed at T = T c for various values of 
the surface - to - bulk coupling ratio Jx/J. Particular at- 
tention will be paid to the crossover regime, which is best 
described by the shape crossover of order parameter and 
energy density profiles. Outside the crossover regime a 
nonasymptotic scaling regime is investigated and charac- 
terized by nonuniversal J\j J dependent exponents. This 
work is mainly focussed on scaling properties and there- 
fore the simulations are restricted to simple cubic lattices 
with (100) surfaces. 



II. MODEL AND SIMULATION METHOD 

The model Hamiltonian describes an isotropic Heisen- 
berg ferro - or antiferromagnet on a simple cubic lattice 
with two (100) surfaces in a cubic geometry. It is given 
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by (see also Ref.i) 

h = -j s >- s i- j i Yl S '- S J' ( 2 - x ) 

where V denotes the set of all lattice sites (volume) and 
A denotes the set of surface sites. The coupling constants 
J and J\ are assumed to have the same sign and the spins 
Si = (Sf,Sf, Sf) with | Si | = 1 only interact with nearest 
neighbors. A nearest neigbor pair of spins is part 
of the interior V\A of the system, if at least one of the 
spins is not part of the surface A. A nearest neighbor- 
pair (i, j) is part of the surface if both spins belong to A. 
The surfaces are defined by the lattice planes (x, y, z = 1) 
and (x, y, z — L), where 1 < x, y < L' . Along the x and 
y directions periodic boundary conditions are applied. In 
order to avoid unwelcomed finite - size effects the geom- 
etry of the lattice is chosen as cubic, i.e., L' = L. The 
cubic shape turns out to be a reasonable compromise be- 
tween achievable system sizes L and the sensitivity of the 
system to the surface - to - bulk coupling ratio J\/ J ■ 

The Monte - Carlo algorithm is chosen as a hybrid 
scheme which consists of Metropolis sweeps, Wolff single 
cluster updatesEa, and overrelaxation sweepstH. Typi- 
cally, one hybrid Monte - Carlo step consists of 10 indi- 
vidual steps each of which can be one of the updates listed 
above. The Metropolis and the Wolff algorithm work the 
standard way, where the reduced coordination number 
of the lattice at the surfaces and the modified surface 
coupling Ji must be taken into account. The acceptance 
probability p of a proposed spin flip in the Metropolis al- 
gorithm is defined by p(AE) = l/[exp(A^/fc B r) + 1], 
where ks is the Boltzmann constant and AE is the 
change in configurational energy of the proposed move. 
The overrelaxation part of the algorithm performs a mi- 
crocanonical update of the configuration by sequentially 
rotating each spin in the lattice such that its energy 
contribution to the energy of the whole configuration 
remains constant. The implementation of this up date 
scheme is straightforward, because according to Eq.(2.1) 
the energy of a spin with respect to its neighborhood has 
the functional form of a scalar product. The angle of ro- 
tation can be chosen randomly for each spin, however, it 
turns out that in view of minimal autocorrelation times a 
reflection, i.e., a rotation of all spins by 180 degrees is the 
most efficient overrelaxation move. One hybrid Monte - 
Carlo step consists of two Metropolis (M) for single clus- 
ter Wolff (C) and four overrelaxation (O) updates. The 
individual updates are mixed automatically in the pro- 
gram to generate the update sequence M O C O C M O 
C O C . The random number generator is the shift reg- 
ister generator R1279 defined by the recursion relation 
X n = X n _ p ®X n _ q for (p,q) = (1279, 1063). Generators 
like these are known to cause systematic errors in com- 
bination with the Wolff algorithmic However, for lags 
(p, q) used here these errors are far smaller than typical 
statistical errors. They Me further reduced by the hybrid 
nature of the algorithm! 4 !]. 



The Monte - Carlo scheme described above is employed 
for lattice sizes L between L = 12 and L = 72. For each 
choice of parameters we perform at least 20 blocks of 10 3 
hybid steps for equilibration followed by 10 4 hybrid steps 
for measurements. Each measurement block yields an es- 
timate for all quantities of interest and from these we ob- 
tain our final estimates and estimates of their statistical 
error following standard procedures. The integrated au- 
tocorrelation time of the hybrid algorithm is determined 
by the autocorrelation function of the energy or, equiva- 
lently, the modulus of the magnetization, which yield the 
slowest modes for the Wolff algorithm. The autocorrela- 
tion times do not exceed 10 hybrid Monte - Carlo steps for 
the largest lattice size (T — T c ), so the equilibration and 
measurement periods quoted above translate to roughly 
100 and 1000 autocorrelation times, respectively. In or- 
der to obtain the best statistics for all magnetic quanti- 
ties a measurement is made after every hybrid Monte - 
Carlo step. All error bars quoted in the following cor- 
respond to one standard deviation. The hybrid scheme 
samples the surfaces of the system reasonably often, so 
a preferential sampling of surface configurations is not 
required. The simulations have been performed on DEC 
alpha and AXP workstations at the Physics department 
of the BUGH Wuppertal. 



III. SURFACE SCALING EXPONENTS 

The simulations presented here have been performed 
at T = T c for several values of the surface - to - bulk cou- 
pling ratip-|Ji/J. The estimate for T c used here is taken 
from Ref.E3, where the critical coupling K c has been de- 
termined as K c = J/{k B T c ) = 0.693035(37). In view of 
the limited system size L < 72 the relative accuracy of 
10~ 4 in T c is sufficient in order to perform a standard 
finite - size scaling analysis in terms of the system size L 
up to the usual corrections to scaling. As our main ref- 
erence for bulk critical exponents we choose the work of 
Guida and Zinn - Justincil, where the bulk critical expo- 
nents of the O(N) universality class have been obtained 
from high order Borcl rcsummcd perturbation theory for 
the Ginzburg - Landau model. 

In order to access surface critical exponents the surface 
magnetization mi, its second moment, and the surface 
energy density e\ has been measured. The layer and 
surface specific heats C\ and C\\ (see Eq.(|l.5|)) have not 
been investgated, because a\ and an can be expressed 
by bulk exponents only. The layer magnetization m(z), 
1 < z < L and the total magnetization m fot are defined 
by 



m(z) = Y S xyz/L 2 
x,y=l 



(3.1) 
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m(z)/L, 



respectively, and m tot = |m tot | denotes the modulus of 
the total magnetization. The magnetization profile m(z) 
is then defined as the projection of m(z) onto the total 
magnetization m tot , i.e., 



m(z) = m(z) • m tot /m tot 



(3.2) 



and mi = (m(l) + m(L))/2 defines the surface magne- 
tization. Note that the surfaces at z = 1 and z = L 
are identical (see Eq.(2.1)). In terms of mtot and m\ the 
layer and surface susceptibilities (i.e., their longitudinal 
components) xi and xn f° r a completely finite system 
are defined as 



Xi = L 2 ((m to tmi) - (mtot){m\))/{k B T) 
Xn =L 2 ((ml) -( mi ) 2 )/(k B T), 



(3.3) 



where (. . .} denotes the thermal average. The energy 
profile e(z) is defined accordingly, where apart from the 
exchange energy between the spins within layer z half 
the interaction energy to the layers z — 1 and z + 1 also 
contributes to e(z), so e tot = X)z=i e ( z )/L i s the total en- 
ergy density. In the following all energies will be giv en i n 
units of ksTc, i.e., extra factors ksT (see, e.g, Eq.flOj)) 
are unity at T = T c . 
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FIG. 1. Surface magnetization m\(L) (see Eq.(3.2)) at 
T — T c as function of the system size L for Ji/J = 0.3 (x, 
solid line), 0.5 (+, short dashed line), 1.0 (*, dashed line), 1.5 
(□, long dashed line), and 2.0 (■, dash-dotted line). Statis- 
tical errors are much smaller than the symbol sizes and the 
various lines are just guides to the eye. For Ji/J < 1 mi dis- 
plays the expected behavior for the ordinary surface univer- 
sality class. For Ji/J = 1.5 the system undergoes a crossover 
towards ordinary surface behavior, whereas for Ji/J = 2.0 
the behavior of mi is inherently nonasymptotic within the 
available range of system sizes. 

The surface magnetization mi at T — T c as function 
of the system size L is shown in Fig.|] for Ji/J = 0.3, 



0.5, 1.0, 1.5, and 2.0. For Ji/J < 1 the functional form 
of m%(L) = (mi) is accurately captured by 



i 1 (L) = B mi L-^(l + C mi L-^), 



(3.4) 



where B mi is the magnetization amplitude and C mi 
is the amplitude of the leading correction to scaling. 

The associated Wegner exponent uj = 0.78 is taken 

from Ref.EJ. A least square fit of Eq.(3.4) to the data 



for Ji/J < 1 displayed in Fig.[l] yields the estimates 
fc/u = 1.185(6), 1.175(13), and 1.171(7) for J x /J = 
0.3, 0.5, and 1.0, respectively. The error indicated 
in parenthesis corresponds to one standard deviation. 
From these estimates one obtaines the weighted average 
0x/v = 1.179(6), where the smallest of the individual 
errors is taken as the final error estimate. Two of the 
three individual estimates are included in the error in- 
terval of the fmal estimate. From the literature value 
v = 0.7073(35)cil one obtains the estimate 



Pi = 0.834(6) 



(3.5) 



for the surface exponent of the magnetization. For 
Ji/J — 1.5 Eq.(3.4) does not capture the functional 
form of mi(L), because within the available range of lat- 
tice sizes the system undergoes a crossover towards the 
asymptotic ordinary surface critical behavior. A more 
detailed discussion of this crossover is postponed to Sec. 
5, where the order parameter and energy density profiles 
are presented. If the decay of mi(L) for J x /J= 1.5 i s 
described by an effective exponent according to Eq. ( |3.4| ) , 
one finds a value around 0.6 for L < 20 and a value 
around 0.9 for L > 48. This indicates that only a part of 
the full crossover process is captured by the simulation. 
This leads to the conclusion that the data for Ji/J = 2.0 
have not yet even entered the crossover regime to ordi- 
nary surface critical behavior. It is insructive to compare 
these data with corresponding data for the Ising model. 
The surface - to - bulk coupling ratio Ji/J — 2.0 already 
belongs to the extraordinary regime of the Ising modelQ, 
where the surface exhibits long - range order at T = T c . 
The comparison is shown in FigJ 2 where the data for 
an Ising model according to Eq. ( |2T[ ) have been obtained 
from a hybrid algortihm which corresponds to the one 
described above, except that c*errelaxation moves can- 
not be performed in this caseEl The surface magne- 
tization decays with an effective exponent of about 0.16 
(see Figj^(a), solid line), whereas for the Ising model (see 
Fig.^(b), solid line) mi(L) approaches the spontaneous 
surface magnetization mio according to 



rai(L) = mio 



(3.6) 



up to corrections to scaling, where fi/v ~ 0.5170 is the 
scaling dimension of the order parameter in the Ising 
universality class and is a nonuniversal amplitude. 
Fig.^ illustrates how the presence of real long - range 
surface order (b) can be distinguished from spurious 
long - range surface order (a) which only appears as 
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a nonasymptotic finite - size effect. However, within a 
typical range of numerically accessible system sizes the 
crossover to the asymptotic ordinary surface behavior 
cannot be observed for J\/J > 2.0 and therefore a data 
analysis within the framework of finite - size scaling only 
yields an effective exponent @i,eff — Pi.effi.Ji/ J)- For 
Ji/J = 2.0 one has /3i, e //(2.0)/V ~ 0.16 (see Fig.|(a)) 
and for J x /J= 3.0 an effective exponent /?i, e / / (3.0) ~ 
0.08 is found (not shown). We do not quote error bars 
here, because the estimates for /3i je //(Ji/J) are pre- 
sumably affected by systematic errors which are larger 
than the statistical ones. In the interior of the sys- 
tem the 'bulk' magnetization = (m(L/2)) obeys 
standard critical finite size scali ng with the exponent 
P/v ~ 0.518EI (see Eqs.Q and Q and Fig.^in Sec. 
5). 
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L 

FIG. 2. Surface magnetization mi(L) for Ji/J = 2.0 at 
T = T c (x) (a) for the Heisenberg model given by Eq.(Er-lt) 
(see also Fig.nl) and (b) for a corresponding Ising modeEI. 
For the Heisenberg model (a) mi(L) decays according to an 
effective power law (solid line) , whereas for the Ising model (b) 
mi(I) increases towards the val ue o f the spontaneous surface 
magnetization according to Eq.(3.6) (solid line). 



The layer susceptibility xi is analyzed in the same way 
as the surface magnetization mi(L). The data are dis- 
played in Fig||. For J\/J<1 the data can be interpreted 
according to 



Xi(L) 



B X1 L^(1 + C X1 L- 



(3.7) 



u> ~ 0.78 as above a least square fit of Eq.(3.7) to the data 
shown in FigJH yields the estimates 71/^ = 1.314(23), 
1.305(12), and 1.308(25) for Ji/J = 0.3, 0.5, and 1.0, 
respectively. As before we adopt the weighted average 
71/f = 1.307(12) as our final estimate, where the small- 
est of the individual errors is taken as the error estimate. 
From the literature value v — 0.7073(35) the estimate 



71 = 0.924(10) 



(3.8) 
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which is the exact analog of Eq.(3.4). The term in paren- 
thesis captures the leading correction to scaling, ji/v is 
the corresponding surface exponent for finite - size scal- 
ing, and B Xl is a nonuniversal overall amplitude. With From Eq.(0) one furthermore obtains the correlation 



FIG. 3. Layer susceptibility xi (see Eq.(gJ)) at T = T c 
as function of the system size L for J\/J = 0.3 (x, solid 
line), 0.5 (+, short dashed line), 1.0 (*, dashed line), 1.5 (□, 
long dashed line), and 2.0 ■ dash-dotted line). Statistical 
errors are smaller than the symbol sizes and the various lines 
are just guides to the eye. For Ji/J < 1 xi displays the 
expected behavior for the ordinary surface universality class. 
For Ji/J = 1.5 the system undergoes a crossover towards or- 
dinary surface behavior, whereas for Ji/J = 2.0 the behavior 
of xi is inherently nonasymptotic within the available range 
of system sizes (see also Fig.hj) . 

is found. For Ji/J > 1.5 the data qualitatively be- 
have as in Figj^: For Ji/J = 1.5 the system undergoes a 
crossover to ordinary surface behavior and for J\j J = 2.0 
the asymptotic behavior is out of reach for the simula- 
tion. Nonetheless, the increase of xi m t ne latter case 
can be described by an effective exponent 7i, e //(</i/>/ = 
2.0)/v ~ 2.31. Likewise, 7i, e //(Ji/J = 3.0) /i/ ~ 2.44 is 
obtained (not shown). Error bars are not quoted for the 
reasons indicated above. 

The exponets [3\ and 71 are no t in dependent. From 
the scaling relations given by Eq.(1.3) and bulk scaling 
relations one can in£er the simple rule {3\ + 71 = (3 + 7. 
From the lit eraturecll one obtains /3 + 7 = 1.7557(56) 
and Eqs.dU) and (|J) yield ft + 71 = 1.758(11) which 
verifies the above scaling law. The effective exponents 
/?i, e// (2.0) and 7i, e //(2.0) for J x /J = 2.0 (see Figs .| and 
g) yield / 9 1)6 / / (2.0) +7i, e //(2.0) ~ 1.75 which is remark- 
ably close to the value of j3 + 7 quoted above. Likewise, 
/3 hef f(3.0) +7i, e //(3.0) ~ 1.78 is found for J x /J = 3.0. 
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exponents 

?7H = 1.358(12) , r}± = 0.697(6) (3.9) 

and the exponent 

7ii = -0.253(9) (3.10) 

of the surface susceptibility xil- At T — T c the surface 
susceptibility therefore behaves according to 



Xii(L) = xno - B X11 L 1 -''^ (1 + C X11 L- 



(3.11) 
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FIG. 4. Surface susceptibility xn at T = T c (x) as function 
of the system size L for (a) Ji/J = 0.3, (b) Ji/J = 0.5, and 
(c) Ji/J = 1.0. Error bars are much sm aller than the symbol 
sizes. The solid lines display fits of Eq.( 3.11 ) to the data (a), 
(b), and (c), respectively. 



TABLE I. Comparison of estimates for /3i from ex- 
periments on Ni(100)li3 (average of three measurements}, 
early Montp-i- Carlo workEj, high temperatipc seriescj, 
field theorjiill, transfer matrix Monte - Carlo£3, massive 
field theoryEl ([0/2], [2/0], [1/1] Pade approximants) , and 
Eq.\"~ 



RefM 


Ref 


30 


Ref. 


H Ref. 


17 


Ref£3 Ref. 


12 


Eq.pl) 


Pi 0.81 


0.75(10) 0.81(4) 0.84(1) 0.80(3) 


0.834(6) 



where the term in parenthesis ca pture s the leading cor- 
rection to scaling. Fits of Eq.( B.ll ) to the data for 
Ji/J< 1 are shown in Fig.[|, where 1 — ?7|| = —0.358 and 



uj = 0.78 are kept fixed and the amplitudes xno, B 



and C Y11 are used as fit parameters. For Ji/J > 1.5 
Eq.( |3.1l[) does no longer describe Xn due to the strong 
nonasymptotic effects described above (for more details 
see Sec. 4). A comparison of different estimates for (3i 
from various sources is shown in Table Q. The estimates 
are rather consistent within their mutual errors, but the 
best agreement is found with the field theoretic estimate 
given in Ref £3. 

We close this section with an investigation of the finite- 
size behavior of the surface energy density e\. The 
data for the dimensionless surface energy density s\ = 
ei/{k B T c ) are shown in Fig.| for 0.5 < Ji/J < 2.0. The 
L dependence of ei for Ji/J = 0.3 (not shown), 0.5 (a), 
and 1.0 (b) is qualitatively different from the L depen- 
dence for Ji/J = 1.5 (c), 2.0 (d), and 3.0 (not shown). 
The leading L dependence of ei is written as 



£x{L) = eio 



(3.12) 
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FIG. 5. Dimensionless surface energy density 
£i = ei/(fcsr c ) at T = T c (x) as function of the system 
size L for (a) Ji/J = 0.5, (b) Ji/J = 1.0, (c) Ji/J = 1.5, 
and (d) Ji/J = 2.0. In (c) and (d) the error bars are much 
sma ller t han the symbol sizes. The solid lines display fits of 
Eq.(3.12) to the data (a), (b), (c), and (d), respectively. 



(2) (3) 

where the amplitudes £io, , and B^ are taken as 
fit parameters. For J±/J < 1 the coefficient Bl t turns 
out to be two orders of magnitude smaller than the coef- 
ficient Bif* which suggests that L~ 3 is the leading finite 
- size correction in this case. According to naive finite - 
size scaling one expects x ei = (1 — ai)/v as the leading 
finite - size exponent for the surface energy density. In 
fact, due to the-s&aling relation ai = a — 1 for the ordi- 



nary transitionliSlia one finds 



(2 - a)/p = d = 3 
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from hyperscaling in agreement with the simulation data 
in Figs.|| (a) and (b). In contrast, B^ and Be? are of 
roughly equal magnitude (and of opposite sign) for the 
case Ji/J > 1.5 displayed in Figs.g (c) and (d). There- 
fore, L~ 2 rather than L~ 3 is the leading finite - size cor- 
rection here. The discrepancy between Figs.[| (a), (b) 
and Figs.[5] (c), (d) is due to strong nonasymptotic con- 
tributions to the finite size behavior for Ji/J > 1.5. One 
possible source of the L -2 contribution is the scaling di- 
mension x f g = (2 — a s )/i> = (2 — a — v)jv = d— 1 = 2 
of the surfac e free energy f s (surface tension) which pre- 
vails in Eq.(3.12) due to the-jnonscaling dependence of 
the surface tension on Ji/Jcj. Another source of the 
leading L~ 2 dependence is the regular (noncritical) finite 
- size behaviorj-Qf the surface tension for periodic bound- 
ary conditionstj. Other finite - size corrections such as 
L^ 1 substantially reduce the goodness - of - fit and can 
therefore be ruled out as leading terms. .Eossible log- 
arithmic corrections of the form L~ 2 logZc3 cannot be 
identified unambigiously. 
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FIG. 6. Transverse component of the surface structure fac- 
tor Sl(p,L) for Ji/J = 0.5 and L = 24 (x), L = 48 (+), 
and L = 72 (*) as function of the momentum tr ansf er p in 
the (100) direction. Numerical evaluations of Eq.(4.2) in the 
limit c — > oo for L = 24, L — 48, and L = 72 (see main text) 
are displayed by the solid, the long dashed and short dashed 
line, respectively. 



IV. SURFACE STRUCTURE FACTOR 

The surface structure factor is given by the discrete 
Fourier transform of the surface spin - spin correlation 
function 

G""(x - x') = {{S" y l S^, yl l ) — (S^ y l ) 2 (4.1) 



where the Fourier transform is taken with respect to 
x— x', x = (x,y), and the upper index a refers to the spin 
component. Note that off - diagonal components of G\ 



vanish identically. The longitudinal and transverse com- 
ponents of G" a with respect to the total magnetization 
are very similar and we therefore restrict the following 
discussion to the transverse component only. The mo- 
mentum transfer p in the (100) direction, i.e., parallel to 
the surface, in units of the inverse lattice constant is given 
by p = 0, 0, 0) = (2irn/L, 0, 0), n = 0, 1, . . . , L/2. Nu- 
merical data of the surface structure factor for Ji/J = 0.5 
in the (100) direction for three different lattice sizes are 
displayed in Fig.||. The data essentially collapse onto a 
single curve for p > 0, finite - size effects arc only visi- 
ble at p — 0, where the structure factor shown in Fig.^ 
reduces to the transvers e com ponent of the surface sus- 
ceptibility xii (see Eq.fl3.llD). For Ji/J = 0.3 identi- 
cal properties are obtained (not shown). It turns out 
that lattice effects in the structure factor near the Bril- 
louin zone boundary can be captured by.-the replacement 
p — > 2 sin(p/2) to a remarkable accuracyC-X The shape of 
the surface structure factor is rea sonab ly well captured 
by the mean - field type expressioncjej 



S?(p,L)= X ? 10 -B° 



2 sin(p/2) 



+ 2sin(p/2) 
2r Q sin(p/2) 



(4.2) 



tanh(2r« sin(p/2)L) 

where the upper index a indicates the component (longi- 
tudinal / or transverse t with respect to m tot ) of 5f (p, L). 
The coefficients Y?m and B" are the coefficients of the 
surface susceptibility Xii (f° r a = I see Eq.Q3.llf))- The 
width parameter r Q and the surface enhancement pa- 
rameter c a are used to fit the momentum dependence 
of the surface stru cture factor and the exponent r)u is 
taken from Eq.( |3.9p . The fit is performed in two stages. 
First, x" ln an d B" are determined from a least square 



fit of Eq.(4.2) for p = to the surface susceptibility. 
Second, the width parameter T a and, if needed (see be- 
low), the surface enhancement parameter c a are adjusted 
to obtain a least square fit of the momentum depen- 
dence to the data. In practice, this fit procedure has 
only been performed for the largest lattice size L = 72. 
For smaller systems the fit parameters are taken from 
L = 72 in or der to test the accuracy of the shape pre- 
dicted by Eq.(4.2). The result for J x /J= 0.5 is displayed 
in Fig^ for the transverse component S{ (p, L) of the sur- 
face structure factor. For the largest system L = 72(*) 
the short dashed line shows the fit, whereas for other 
lattice size s su ch as L — 24(x) and 48(+) the evalua- 
tion of Eq.flOl) for L — 24 (solid line) and L = 48 (long 
dashed line) confirms the predicted shape to a remark- 
able accuracy. For J\/J — 0.3 and 0.5 the parameter 
c a is very large and therefore the surface enhancement 
prefactor (c a — 2 sin(p/2))/(c Q + 2 sin(p/2)) can be omit- 
ted. Noticeable deviations between the data and Eq. fl4~2]) 
for a = t only occur for small p on the sm alle r lattices, 
where additional lattice corrections to Eq.( |4.2j ) may be- 
come important. However, for any p > the data for all 
system sizes collapse onto a single curve, which can be 
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obtained from Eq. ( |4.2| ) by performing the limits c" — > oo 
and L — > oo at finite p. 

For Ji/J= 1.0 the surface enhancement parameter c a 
becomes important. The data are shown in Fig.^, where 
the surface enhancement prefactor provides an important 
correction to the momentum dependence of S1{p,L). 
The fit procedure works as described above, but devi- 
ations from the assumed shape for small p also occur for 
L = 72. The surface structure factor is more sensitive to 
crossover phenomena occurring for larger values of J\j J 
than the surface quantities discussed in the previous sec- 
tion. Corrections to scaling, which are not included in 



Eq.(4.2), may also account for part of the deviations be- 
tween the data and the simple model for the shape func- 
tion. Note that the analytic results obtained in Refr3 do 
not hold for the cubic geometry used here. Furthermore, 
nonasymptotic surface enhancement corrections, which 
have not been considered in Ref.E3, become essential for 
the data analysis. For p > the data still collapse onto 
a single curve, which is given by Eq.(4.2) in the limit 
L — > oo at finite p and finite c a ~ 15. 
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FIG. 7. Transverse component of the surface structure fac- 
tor S{(p,L) for Ji/J = 1.0 and L = 24 (x), L = 48 (+), and 
L = 72 (*) as function of the momentum t rans fer p in the 
(100) direction. Numerical evaluations of Eq.( |4.2| ) for L — 24, 
L = 48, and L = 72 (see main text) are displayed by the solid, 
the long dashed and short dashed line, respectively. Surface 
enhancement corrections are essential for the data analysis. 

The behavior of S[(p,L) for Ji/J = 2.0 is shown in 
Fig.|8|, for Ji/J — 3.0 similar results have been obtained 
(not shown). Although scaling appears to be valid to 
a very high degree of accuracy, the behavior of S{(p,L) 
is very different from Figs.^ and [t]. First, the surface 
susceptibility grows according to the power law 



-n\\,eff(Jl/J) 



(4.3) 



rather than approaching a fi nite limit as in Eq.(3.11). 
From a least square fit of Eq.(4.3) to x\\ one obtains the 
effective exponent f]\\,eff(Ji/J = 2.0) ~ —0.64, where 
deviations from the pure power law given by Eq. (|4.3|) are 



very small. For J\/J — 3.0 one obtains r)\\ te f f{%.Q) ~ 
—0.82. We refrain from quoting error bars here, because 
the values for V\\,eff(Ji/ J) mav be affected by systematic 
errors due to corrections to Eq.(4^2) of unknown form. 
Mutual interaction between the surfaces mediated by the 
bulk may also cause systematic errors. With B\ w and 
f]\\,eff{Jl/J) taken from Eq.( [4.3| ) the shape function 



-2sin(p/2) 
X11 c« + 2sin(p/2) 

2r Q sin(p/2) 



(4.4) 



"1 ')ii,e//("'i/- 7 )- 1 



tanh(2r Q sin(p/2)L) 
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FIG. 8. Transverse component of the surface structure fac- 
tor Sl{p,L) for Ji/J = 2.0 and L = 24 (x), L = 48 (+), 
and L = 72 (*) as function of the momentum tran sfer p in 
the (100) direction. Numerical evaluations of Eq.(4.4) for 
L = 24, L = 48, and L — 72 (see main text) are displayed by 
the solid, the long dashed and short dashed line, respectively. 
The behavior of Si is is inherently nonasymptotic. 

is fitted to S\, where the remaining two parameters 
c a and r a are used. As described above the fit is only 
performed for L — 72 (short das hed line in Fig.||). For 
L = 24(x) and L = 48(+) Eq.Q (a = t) is shown 
for L = 24 (solid line) and L = 48 (long dashed line), 
where all fit parameters are kept fixed. The shape func- 
tion given by Eq.(4.4) is remarkably accurate. How- 



ever, the observed scaling is completely different from 
the asymptotic scaling shown in Fig.^| and, apart from 
surface enhancement corrections, in FigJ^. The effective 
exponent r}\\ e ff(J\/ J) is not related to the surface ex- 
ponent 77|| given by Eq.(3.9), because it is nonuniver- 
sal, i.e., it depends on J\/J. From general considera- 
tions it is tempting to pose the (effective) surface scal- 
ing law /3i, e //(Ji/J) = v\d - 2 + V\\,eff {Ji/J)]/2 (see 
Eq.(1.3)). The direct determination of Pi, e ff from m\{L) 
for Ji/J> 2.0 is plagued with considerable uncertainties, 
because the decay of m\(L) with L becomes quite slow 
for Ji/J > 2.0. Systematic errors may exceed the for- 
mal statistical error of a fit in this case and therefore 
the aforementioned effective scaling law cannot be con- 







firmed unambigiously. Note that the crossover regime to 
the asymptotic scaling remains unaccessible within the 
range of system sizes used here. Nonetheless, the accu- 
racy of the scaling law for Si (p, L) for the nonasymptotic 
regime still lacks theoretical understanding. 

For Ji/J = 1.5 none of the above descriptions applies 
to the data and scaling appears to be violated. This is 
in accordance with the findings of Sec. 3, where m\{L) 
and Xl(L) display sizable deviations from simple scaling 
laws. The effects of the crossover to the ordinary sur- 
face universality class are particularly pronounced in the 
shape function of the order parameter profile to which 
we turn in the following section. 



V. PROFILES 



The order parameter profile m(z) (see Eq.(3.2)) pro- 
vides local information about the order in the system. 
Furthermore, {m(z)} is by construction very sensitive to 
the system size at T = T c and should therefore be a 
valuable probe for scaling behavior which is easier to in- 
terpret than the structure factor. For Ji/J < 1 all previ- 
ous investigations have shown that the system essentially 
displays the asymptotic critical behavior of the ordinary 
surface universality class. The magnetization profile con- 
firms this again, so the scaling analysis can be restricted 
to the case J\/ J — 0.5. The scaled magnetization profile 
M(z/L) is defined by 



M(z/L) = {m{z))/m b , 



(5.1) 



where 



(m(L/2)) = B mi L- ' u (1 + C mb L-") . (5.2) 



The factor in parenthesis captures corrections to scaling 
and the exponents 3, v, and u are taken from Rcf.Eil. 
The co effic ient C mb is determined by a least square fit 
of Eq.(5.2) to mi,. For numerical convenience z — 1/2, 
z = 1, . . . , L is chosen as the position coordinate for the 
profile. The scaling plot for J x /J= 0.5 is shown in Fig||, 
where z now refers to the shifted layer index. The data 
scale very accurately and the shape of the profile is in 
accordance with the expectation for the ordinary tran- 
sition. The scaling function M(z/L) can be represented 
by the simple fit formula 



M(C)=B M [( C + Co )(l_ C + Co)] 



(ft-/?)/" 



(5.3) 



where Co = zq/L and zp-is the extrapolation length, 8 an d 
v are taken from Ref.cil, and B\ is taken from Eq.(3.5). 
From a least square fit of Eq^x^) to the data for L = 72 
one finds z ~ —0.26 in units of the lattice constant. For 
Ji/J = 0.3 and 1.0 one finds z — —0.34 and z ~ 0.46, 
respectively. The choice 0.3 < Ji/J < 1 therefore yields 
quite accurate realizations of Dirichlet boundary condi- 
tions, which characterize the ordinary surface universal- 
ity class. 
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FIG. 9. S calin g function M(z/L) of the magnetization pro- 
file (see Eq.Q)) at T = T c for J t /J = 0.5 and L = 24 (x, 
solid line), 36 (+, short dashed line), 48 (*, dashed line), 60 
(□, long dashed line), and 72 (■, dash-dotted line). Statis- 
tical errors are much smaller than the symbol sizes and the 
various lines are just guides to the eye. Scaling is obtained in 
accordance with the ordinary surface universality class. 
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FIG. 10. Scaling function M(z/L) of the magnetization 
profile (see Eq. ([□])) at T = T c for Ji/J = 2.0 and L = 24 (x , 

dashed line), 60 
Statistical 



solid line), 36 (+, short dashed line), 48 (t 
(□, long dashed line), and 72 (■, dash-dotted line), 
errors are much smaller than the symbol sizes and the various 
lines are just guides to the eye. The surface magnetization is 
strongly enhanced compared to the bulk magnetization (spu- 
rious long - range surface order). 



For Ji/J = 2.0 the application of Eqs.(M) and ( |5.2j ) 
to the data yields the result shown in Fig.Hj Except at 
the surface layers scaling is fulfilled very accurately, how- 
ever, the shape of m(z) does not show the expected fixed 
point form. Instead, a strong enhancement of the surface 
magnetization over the magnetization in the interior is 
obtained. A behavior like this is typical for the extraor- 
dinary transition which does not occur for the Heisenberg 
model defined by Eq.(2.1) in d — 3, i.e., the system dis- 
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plays spurious long-range surface order. If the lattice 
size L could be increased further one should therefore 
find a crossover to the ordinary behavior displayed in 
Fig.|| (see Fig.||). Nonetheless, the shap e of the scaling 
function M(£) is again captured by Eq.(5.3), if the ex- 
ponent is used as a third fit parameter. It is tempting 
to write this exponent in the form (f3i, e ff(Ji/J) — 
however, the result for f3\ je f f{J\/ J) obtained this way is 
not compatible with the estimate obtained from mi(L), 
which may be due to systematic errors of various kinds 
discussed above. 
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FIG. 



1 1. S caling plot of the magnetization profile M(z/L) 
(see Eq.(5.1)) at T = T c as function of z/L for J\/J = 1.5 
and L = 24 (x, solid line), 36 (+, short dashed line), 48 (*, 
dashed line), 60 (□, long dashed line), and 72 ■ dash-dotted 
line). Statistical errors are much smaller than the symbol 
sizes and the various lines are just guides to the eye. Scaling is 
grossly violated durin g th e crossover from the state of spurious 
surface order [see FigJlOJ) to the ordinary surface universality 
class (see Fig.^) as L is increased. 

The question how the system actually performs the 
crossover is answered in Fig. 11, where M(z/L) is shown 
for Ji/J — 1.5. For smaller systems the surface magneti- 
zation is still enhanced over m^, but as L is increased, the 
maximum of the profile at z = L/2 finally exceeds mi and 
the profile shape approaches the fixed point form shown 
in Fig.|^. The influence of the surface coupling on (m(z)) 
is confined to the two outermost lattice layers on either 
side of the cube, whereas the curvature of the remainder 
of the profile already has the "correct" sign for all lattice 
sizes shown in Figjll]. For L < 24 the profile becomes 
flatter in the middle and the profile shape approaches the 
one displayed in Fig. 10. Note that ev en a t z = L/ 2 the 
data fail to collapse according to Eqs.(5.1) and fl5.2p. Al- 



though the isotropic Heisenberg model does not display 
long - range surface order in the thermodynamic limit, 
on a finite lattice a strong enhancement of the surface 
magnetization mi over the bulk magnetization mj does 
occur for sufficiently strong surface couplings as a finite 
- size effect for a certain range of system sizes. 
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z/L 

FIG. 12. Sc alin g function E(z/L) of the energy density 
profile (see Eq.(Q) in units of k B T c at T = T c for J x /J = 0.5 
and L = 24 (x, solid line), 36 (+, short dashed line), 48 (*, 
dashed line), 60 (□, long dashed line), and 72 (■, dash-dotted 
line). Statistical errors are smaller than the symbol sizes and 
the various lines are just guides to the eye. Scaling is obtained 
in accordance with the ordinary surface universality class. 

The scaling properties of the energy density profile are 
a little more delicate, because a background energy den- 
sity must be subtracted from the profile in order to obtain 
scaling. One finds the scaling form 



E(z/L) = ((e(z))-e )/L- 



-{l-ot)/v 



(5.4) 



where E{C/) is the scaling function and corrections to 
scaling have been disregarded. For numerical conve- 
nience z — 1/2, z = 1, . . . , L is again chosen as the po- 
sition coordinate for the profile. The scaling plots for 
0.3 < Ji/J < 1.0 are well represented by the scaling plot 
for Ji/J = 0.5 which shown in Fig.[l2[ As in Figsj9| - 
[n] z refers to the shifted layer index. The shape of the 
scaling function is as expected for the ordinary univer- 
sality classed. For L > 24 the data collapse rea sonably 
well which confirms scaling according ia Eq.(5.4), where 
(1 — ot)/v ~ 1.586 according to RefJlil. The shape of 
the scaling function EJX) can be approximated by the fit 
formula (see also Ref.c3) 



E(0 = Be 



it/ sin 



. C + Co 

l + 2Co 



1 



(5.5) 



where Co = z o/L is the scaled extrapolation length of 
the profile (see Eq.( |5.3D ). For small arguments £ = z/L 
Eq.( |5.5|) captures the algebraic increase of (e(z)) cor- 
rectly and the extrapolation length zq beco mes negligibly 
small. However, away from the surface Eq.(5.5) captures 
the shape function of the energy density profile only in a 
qualitative sense. 
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E(z/L) 




z/L 

FIG. 13. Scaling function E(z/L) of the energy density 
profile (see Eq.(Q) in units of k B T c at T = T c for J x /J = 2.0 
and L = 24 (x, solid line), 36 (+, short dashed line), 48 (*, 
dashed line), 60 (□, long dashed line), and 72 (■, dash-dotted 
line). Statistical errors are smaller than the symbol sizes and 
the various lines are just guides to the eye. The surface en- 
ergyis strongly enhanced compared to the bulk energy (see 
Fig©. 



extrapolation length zq ~ —0.12 (in units of the lattice 
constant) is still very small. The behavior the energy 
density profile displayed in Figjl3| is typical for a system 
at the extrordinary transition which does not exist for 
the Heisenberg model in d = 3 with nearest neighbor in- 
teractions. According to the above discussion the shape 
of the energy density profile given by Fig.[l3| is governed 
by the presence of spurious long - range order in the sur- 
face. The crossover to the aymptotic shape (see Figjl^) 
will occur if L is increased further. For Jx/J = 2.0 the 
crossover regime is out of reach, but for Jx/J = 1.5 this 
crossover takes place within the range of accessible sys- 
tem sizes as shown in Fig.|l4|. As in Fig.[ll] scaling is vio- 
lated. The curvature of the profile near z = L/2 changes 
sign between L — 36 and L — 60 and for L > 60 the 
profile approaches its asymptotic shape. The nonasymp- 
totic surface effects are more pronounced here than for 
the magnetization profile and penetrate deeper into the 
system, but the magnitude of the surface induced en- 
hancement of the energy density at the surface decays 
quickly with increasing L. The same crossover behavior 
can be observed for the critical Ising model slightly below 
the SB multicritical point for, e.g., Jx/J = 1.45e3. 



VI. SUMMARY AND OUTLOOK 



E(z/L) 




z/L 

FIG. 1 4. S caling plot of the energy density profile E(z/L) 
(see Eq.(5.4)) in units of ksT c at T = T c as function of z/L 
for Ji/J = 1.5 and L = 24 (x, solid line), 36 (+, short 
dashed line), 48 (*, dashed line), 60 (□, long dashed line), 
and 72 (■, dash-dotted line). Statistical errors (not shown) 
are slightly larger than the symbol sizes and the various lines 
are just guides to the eye. Scaling is grossly violated during 
the crossover from the state of enhanced surface energy (see 
Fig.[l3|) to the ordinary surface universality class (see Fig.|l^) 
as L is increased. 

For Jx/J = 2.0 scaling of the data according to 
Eq.( |5.4| ) works even better as shown in Fig.[ll| Note 
that the value of the reference energy density e does not 
depend on Jx/J- The data for system sizes L > 24 col- 
lapse onto a single curve, which is represented by Eq.(5.5) 
with a much better accuracy than for Jx/J < 1.0. The 



In the absence of symmetry breaking fields the asymp- 
totic critical scaling behavior of surfaces of a critical 
d = 3 dimensional Heisenberg magnet with short - range 
interactions is always governed by the ordinary surface 
universality class. For finite systems, however, a cru- 
cial interplay between the available system size and the 
value of the surface - to - bulk coupling ratio Jx/J de- 
termines whether or not the asymptotic surface scaling 
behavior can actually be observed. Within the range of 
system sizes 12 < L < 72 for the L x L x L geome- 
try used here critical behavior in the ordinary surface 
universality class can be observed for Jx/J < 1-0. The 
scaling exponents found numerically in this case are con- 
sistent with rigorous scaling laws and estimates for the- 
ses exponents found previously by various analytical and 
numerical methods. The shape of the surface structure 
factor Si (p, L) is captured by very simple mean - field 
like expressions, in which two amplitudes are fixed by 
a fit to the surface susceptibility for p ~ 0. A width 
parameter and, if needed, a surface enhancement param- 
eter then determine the shape of the momentum depen- 
dence. It turns out, that finite - size and lattice effects 
are very accurately described by the pseudo scaling ar- 
gument 2 sin(p/2)L which replaces the true scaling argu- 
ment pL. These properties of Sx{p,L) are essential for 
the data interpretation of the dynamic surface structure 
factor which is the key quantity for the interpretaion of 
neutron scattering data on magnetic surfaces and will 
therefore be the main focus of ensuing work. The order 
parameter and energy density profiles are less relevant 
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for experiments, but they are easier to interpret and yield 
valuable insight into the scaling behavior of the system. 
Either profile is found to scale in accordance with the 
ordinary surface universality class. 

For Ji/J > 2.0 scaling is still found for all quanti- 
ties under investigation, however, the scaling exponents 
are replaced by effective ones and their values depend on 
Ji/J. The scaling relations cannot be verified unambi- 
giously, because some of the numerical estimates for the 
effective exponents are presumably affected by system- 
atic errors of unknown magnitude. Such errors may ensue 
due to unknown corrections to the effective scaling laws 
or due to an effective interaction between the two surfaces 
mediated by the bulk system in between. Nontheless, the 
effective scaling properties of the surface structure factor 
S\ (p, L) provide valuable information for the analysis of 
its dynamic counterpart. The striking scaling properties 
found here still await theoretical explanation. 

Coupling ratios J\/J > 2.0 are too large to access the 
crossover regime from the state of enhanced surface mag- 
netization (spurious long - range surface order) to the 
asymptotic (ordinary) surface scaling. For Ji/J = 1.5, 
however, this crossover becomes the dominating feature 
in the finite - size behavior of all quantities under investi- 
gation, at least within the range of system sizes used here. 
The value Ji/J = 1.5 only marks the location of the 
crossover regime for the system sizes at hand rather than 
a sharp transition in the surface behavior of the Heisen- 
berg model. In the crossover regime scaling is violated 
and further theoretical insight is needed for a purposeful 
data analysis. Scaling of bulk quantities is also affected 
by particularly large correction terms. The crossover pro- 
cess itself is best visualized in the shape crossover of the 
magnetization and the energy density profiles, which oc- 
curs as L is increased for J 1 /J= 1.5 and T = T C . 
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